% GAUTIER LE BIHAN - 2020
% Replication files for "Shocks vs Menu Costs: Patterns of Price Rigidity in an Estimated Multi-Sector Menu-Cost Model?" Review of Economics and Statistics
%
% This code produces Table 3
clear;
tic
cd ..\..\Simulations_VC\MS_produits_MC

load ..\..\Simulations_VC\MS_produits_MC\stat4b
param=stat4b;

load ..\..\Estimation_param\MS_produits_MC\actual_moments_k
weight=actual_moments_k(:,12);
sum(weight)

secteur=actual_moments_k(:,1);

p0=param(:,1);
%muc=param(:,2);
sig_ea=param(:,3);
rho_a=param(:,4);
sig_a=sig_ea./sqrt(1-rho_a.^2);
share_calvo=param(:,end-3);
muc=param(:,end-2);
av_muc=param(:,end-1);

sig_tilda=sqrt((2*sig_ea.^2)/(1+0.6946))
mu_tilda=muc.*(1-share_calvo);
stat=[p0 muc av_muc sig_ea rho_a sig_a share_calvo sig_tilda mu_tilda];
%mean(stat)
%median(stat)
size_param=size(stat,2)

stat_q1=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_q1(i) = mean(sortx([j j+1]));
            else
                
                stat_q1(i) = sortx(j+1);
            end
     stat_q1(end)= csumw(end);
     stat_q1(end-1)= v;
    end  
 end


stat_med=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/2;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_med(i) = mean(sortx([j j+1]));
            else
                
                stat_med(i) = sortx(j+1);
            end
     stat_med(end)= csumw(end);
     stat_med(end-1)= v;
    end  
 end
%"stat_med"
%stat_med


stat_q3=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = 3*sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_q3(i) = mean(sortx([j j+1]));
            else
                
                stat_q3(i) = sortx(j+1);
            end
     stat_q3(end)= csumw(end);
     stat_q3(end-1)= v;
    end  
 end



meanw_data_all=weight'*stat/sum(weight);
for i=1:size(stat,2)
stdw_data(i)=sqrt(weight'*((stat(:,i)-meanw_data_all(i)).^2)/sum(weight));
end
fprintf("Stats table Estimates - PANEL C")
[meanw_data_all(1:size_param);
stat_med(1:size_param);
stat_q3(1:size_param)-stat_q1(1:size_param);
stdw_data(1:size_param)]
%'interq / med'
%(stat_q3-stat_q1)./meanw_data_all
%(stat_q3-stat_q1)./stat_med


%IDEM Excluding energy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
stat=stat((stat(:,1)<=0.3),:);
weight=weight((stat(:,1)<=0.3),:);

stat_q1=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_q1(i) = mean(sortx([j j+1]));
            else
                
                stat_q1(i) = sortx(j+1);
            end
     stat_q1(end)= csumw(end);
     stat_q1(end-1)= v;
    end  
 end


stat_med=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/2;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_med(i) = mean(sortx([j j+1]));
            else
                
                stat_med(i) = sortx(j+1);
            end
     stat_med(end)= csumw(end);
     stat_med(end-1)= v;
    end  
 end

stat_q3=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = 3*sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_q3(i) = mean(sortx([j j+1]));
            else
                
                stat_q3(i) = sortx(j+1);
            end
     stat_q3(end)= csumw(end);
     stat_q3(end-1)= v;
    end  
 end


meanw_data_all=weight'*stat/sum(weight);
for i=1:size(stat,2)
stdw_data(i)=sqrt(weight'*((stat(:,i)-meanw_data_all(i)).^2)/sum(weight));
end

fprintf("Idem excluding energy")

[meanw_data_all(1:7);
stat_med(1:7);
stat_q3(1:7)-stat_q1(1:7);
stdw_data(1:7)]
'interq / med'
(stat_q3-stat_q1)./stat_med


